Extra-Large Remnant Recoil Velocities and Spins from 
Near-Extremal-Bowen- York-Spin Black-Hole Binaries 
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We evolve equal-mass, equal-spin black-hole binaries with specific spins of a/iriH ~ 0.925, the 
highest spins simulated thus far and nearly the largest possible for Bowen-York black holes, in a set 
of configurations with the spins counter-aligned and pointing in the orbital plane, which maximizes 
the recoil velocities of the merger remnant, as well as a configuration where the two spins point in the 
same direction as the orbital angular momentum, which maximizes the orbital hang-up effect and 
remnant spin. The coordinate radii of the individual apparent horizons in these cases are very small 
and the simulations require very high central resolutions (h ~ M/320). We find that these highly 
spinning holes reach a maximum recoil velocity of ~ 3300 km s _1 (the largest simulated so far) 
and, for the hangup configuration, a remnant spin of a/mn ~ 0.922. These results are consistent 
with our previous predictions for the maximum recoil velocity of ~ 4000 km s _1 and remnant 
spin; the latter reinforcing the prediction that cosmic censorship is not violated by merging highly- 
spinning black-hole binaries. We also numerically solve the initial data for, and evolve, a single 
maximal-Bowen- York-spin black hole, and confirm that the 3-metric has an 0(r~ 2 ) singularity at 
the puncture, rather than the usual 0(r~ 4 ) singularity seen for non-maximal spins. 

PACS numbers: 04.25.Dm, 04.25. Nx, 04.30.Db, 04.70.Bw 



I. INTRODUCTION 

Highly spinning black holes play an important role 
in astrophysics, from powering active galactic nuclei 
(AGN), to 7-ray bursts (GRB) and quasars. While 
the (indirect) observational evidence for the existence of 
black holes, stellar mass or supermassive, is overwhelm- 
ing, the actual observational evidence for black-hole spin 
is scarce. There have been attempts to measure the cen- 
tral black hole spin in AGN P, H|, Seyfert galaxies 0, 
and quasars The x-ray spectra of accretion disks 
around stellar mass black holes can also provide infor- 
mation about their spins [!, H, 0| • 

The recent dramatic breakthroughs in the numerical 
techniques to evolve black-hole-binary spacetimes [1, H, 
[To| has led to rapid advancements in our understanding of 
black-hole physics. Notable among these advancements 
are developments in mathematical relativity, including 
systems of PDEs and gauge choices fllL lla 13 1, the ex- 
ploration of the cosmic censorship [5ll5llm ll7.1l8l, and 
the app lication of isolated horizon formulae 15l . la , liH 
|20| . 1 2 it [2^ | . These breakthroughs have also influenced 



the development of data analysis techniques through 
the matching of post-Newtonian to fully-numerical wave- 
forms [23l I24I I25I |26| . In particular, the moving punc- 
tures approach proved to work in a wider realm than 
was originally thought. Notably, it has been successfully 
applied to many-black-hole spacetimes W\ l28l. and to 
black-hole-neutron star evolutions [IE l3dl. l3ll. |33| . Sim- 
ilarly, the recent discovery of very large merger recoil 



kicks [33J, |3J, |35j, |36j, |3J, |3J, |3J, |40j has had a great im- 
pact in the astrophysical community, with several groups 
now seeking for observational traces of such hig h speed 
holes as the byproduct of galaxy collisions [HI, E2 EE EH • 

The first study of generic black-hole-binary configura- 
tions (i.e. binaries with unequal component masses and 
spins, and spins not aligned with each other or the or- 
bital angular momentum) was described in Ref. [33j , and, 
based on the results of that study, a semi-empirical for- 
mula to estimate the recoil velocities of the remnant 
black holes was proposed, finding recent confirmation 
in [35|, EE EE| • The spin contributions to the recoil ve- 
locity are generally larger than those due to the unequal 
masses, and, in particular, the spin component in the or- 
bital plane has the largest effect (33[, leading to a maxi- 
mum recoil velocity of about 3500 — 4000 km s _1 [35| . 
The recoil velocities acquired by the remnant of the 
merger of black-hole binaries has many interesting as- 
trophysical consequences (33[, particularly since spin- 
ning black holes can accelerate the merged hole high 
enough to eject the remnant from the host galaxy. Re- 
cently a quasar, displaying blue shifted emission lines by 
2650 km/ s with respect to its host galaxy, has been ob- 
served 

In all the above simulations, the evolution was started 
using conformally flat initial data. This choice has the 
advantage of being easy to implement, with the (ap- 
parently) minor drawback of introducing a short, non- 
physical burst of radiation at the start of the simula- 
tion. Apart from this initial burst, there appears to 
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be no unphysical behavior associated with conformally 
flat initial data, and this choice remains popular (See 
Refs. [HI H§] for more astrophysically realistic initial 
data using post-Newtonian information). The simplest 
initial data, Bowen-York (BY), gives the extrinsic curva- 
ture, Kij analytically by assuming that it is transverse 
and traceless (See (5 Of for an alternative prescription). 
A particularly interesting feature of the conformally flat 
ansatz for the 3-metric is that these data cannot model 
maximally rotating Kerr black holes (such holes are not 
conformally flat in any smooth slice), but have a limit- 
ing spin Ha HU of S/M\ DM fa 0.928 for BY data and 
S/M\ n n/f ~ 0.932 for conformally Kerr extrinsic curva- 
ture [50| • Here S denotes the spin of the black hole and 
Madm the total ADM mass. 

The spin of the merger remnant is similarly important 
both because high-spin black holes are more efficient at 
converting accreting matter into radiation than lower- 
spin holes, and because of open questions regarding cos- 
mic censorship. This issues have already been studied 
in the 'Lazarus' approach to numerical evolutions [52| 
and in early evolutions using the 'moving punctures' ap- 
proach 0, HE Gil ■ Recently the issue has been revisited 
in the context of unequal mass holes [l7l Il8l l53l l54l f55j] . 
and highly-elliptical, equal-mass binaries [56L I57T] . In the 
current work we show that, for the maximum possible 
spin for Bowen-York black holes, the merger remnant will 
always have sub-maximal spins. 

This paper is organized as follows: in Sec. Ill Al we de- 
scribe the Bowen-York initial data for a single spinning 
black hole and how one can obtain the maximum pos- 
sible spin for a Bowen-York black hole, in Sec. Ill CI we 
describe how we obtain initial data for black-hole bina- 
ries with nearly maximal BY spin, in Sec. lIIII we describe 
the numerical techniques used to evolve these binaries, in 
Sec. HVI we give the results and analysis from the numer- 
ical evolutions, and in Sec. fVl present our conclusions. 



II. INITIAL DATA 

An initial data set for the Einstein vacuum equations 
is given by a triple (M , jij , ) where M is a connected 
3-dimensional manifold, 7^ a (positive definite) Rieman- 
nian metric, and a symmetric tensor field on M. , such 
that the constraint equations 



DjK lj - D l K = J\ 
R-KijK^ + K 2 = 2/i, 



(1) 
(2) 



are satisfied on M.. Where D and R are the Levi-Civita 
connection and the Ricci scalar associated with 7^ , and 
K = KijjV. The vector ,P and the scalar function \x 
are determined by the stress energy tensor of the sources 
which describes the matter content of the spacetime. In 
these equations the indices are moved with the metric 7^ 
and its inverse r ) tJ . 



A. Maximum Spin Bowen-York initial data 

The Bowen-York family of initial sets [58| represents 
a relevant class of data suitable for numerical simula- 
tions of black-hole binaries. They are constructed using 
the conformal method for solving the constraint equa- 
tions ((I])-® (for a recent review on this method see [59| 
and references therein). Let us consider one member of 
this family, namely the spinning Bowen-York data. These 
data describe a (non-stationary) black hole with intrin- 
sic angular momentum. In order to construct the data 
we prescribe K %3 , a symmetric, trace free and divergence 
free tensor with respect to the flat metric &y 
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K ij = —n {l e j)kl S k 



ni , 



(3) 



where r is the spherical radius, n % the corresponding ra- 
dial unit normal vector, e l ^ k the flat volume element and 
Sk an arbitrary constant vector which will give the total 
spin of the data. In this equation the indices are moved 
with the flat metric <5.y . The data are given by 



(4) 



where the conformal factor satisfies the following equa- 
tion (which is the conformal version of the Hamiltonian 
constraint ([2|)) 



A</> = - 



18S 2 sin 2 6 
8rV ' 



(5) 



where S 2 = SiSjS 1 ^ and A is the flat Laplacian. 

For any solution (p of Eq. ([5]) the metric 7.^ and 
given by define a solution of the vacuum (i.e. J 1 — 
/i = 0) constraint equations |T])-(E]). in order to find a 
unique and physically relevant solution of ([5]) we need to 
impose appropriate boundary conditions for ip. This is 
essentially equivalent to prescribe the manifold M. of the 
initial data. For example, Eq. ([5|) is singular at r = 0, it 
follows that the solution ip can not be regular at r = and 
hence the origin can not be in the manifold M.. That is 
M 3 is not allowed as manifold in this class of initial data. 
In the present case the manifold will be M. = K 3 \ {0}, 
the origin r = represents another end of the initial data. 

Boundary conditions for black holes are prescribed as 
follows. Let mP > be an arbitrary number. Define the 
function u by 
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m 
~27 



(0) 



Using ([5]) we obtain an equation for u outside the origin 

18rS 2 sin 2 9 



(rlf + ru) 



If u is positive (this will be the case as a consequence of 
the maximum principle) then the denominator of Eq. ([7]) 
never vanishes and hence this equation is regular at the 
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origin. The idea is to impose this equation also at the 
origin, that is, we want to solve ([7]) in R 3 with a boundary 
condition at infinity 



lim u = 0. 



(8) 



It is well known that for each m p > there exists a 
unique solution u of (J7J which satisfies the boundary 
condition ([5]). The solution is positive and from stan- 
dard elliptic theory it follows that it is smooth outside 
the origin and it is C 2 at the origin (this loss of differen- 
tiability is due to the presence of the function r on the 
right hand side of Eq. (UJ which is not smooth at the 
origin). This is what in the numerical relativity is called 
the puncture method [60| , the parameter mP is called the 
mass parameter of the puncture. By Eq. ((S]), the singu- 
lar part of ip at r = is C(l/r), this implies that the 
physical fields 7y , are asymptotically flat at the end 
r = 0. 

The physical parameters of the data are given by S 
which represent the angular momentum of the data and 
the total ADM mass Madm, which is given by the fol- 
lowing formula 



Madm = m p + m u , 
where m u is the coefficient 0(l/r) of u, that is 



(9) 



(10) 



asr-> oo. Note that in order to calculate the total mass 
Madm we need to solve the non-linear equation ([7]) 

The solution u depends on the coordinates x and the 
two parameters m p and S. However, there exists a scale 
invariance for equation ((7J . Namely, if we have a solution 
u(m p ,S,x), for parameters S and m p , then the rescaled 
function u(Xm p , X 2 S, Ace), where A is an arbitrary positive 
number, is also a solution. This means that the solution 
depends non trivially only on one parameter. We chose 
to fix S and vary m p . Note that the following quotient 
is scale invariant 



M ADM 



(11) 



For a Kerr black hole we have £5 < 1 and eg = 1 im- 
plies that the black hole is extreme. For general axi- 
ally symmetric vacuum black holes (which in particu- 
lar includes the spinning Bowen-York data) we also have 
es < 1 61] [62] and £5 = 1 if and only if the data are 
slices of extreme Kerr black hole. Since the Bowen-York 
data are not slices of the extreme Kerr black hole it fol- 
lows that es < 1 for this family. What is the maximum 
value for es in this family? This question was explored 
numerically in [5(| [HHHf. In these references it was ob- 
served that in the limit S — > 00 (for fixed m p ) the ratio 
es reach an asymptotic maximum value. Because of the 
scaled invariance mentioned above this limit is equivalent 
to m p — > 0, with S fixed. What was not clear at all is 



that in fact in the limit we get a well behaved solution of 
the constraint equation. This is precisely the question we 
want to answer here. That is, we want to give numerical 
evidence that the limit 



u(S, x) = lim u(m p , 

mP— *0 



S,x) 



(12) 



exists and defines a solution of the constraint equations. 
An analytical proof of this is detailed in [64| . We will call 
this new solution the extreme Bowen-York spinning data 
because it has the maximum amount of angular momen- 
tum per mass in this family (92j. 

In the limit m p — ► the difference between tp and u 
is just a constant and Eq. {7} become singular at the 
origin, hence the limit solution u can not be regular at 
the origin. If we assume that u = 0(r a ) at r = for some 
real number a, using Eq. ([7J we get that a — —1/2. That 
is, we expect the following behavior at the origin 



O 



1 



(13) 



This behavior is confirmed by the numerical simulations 
presented in the next subsection III Bl 

We have seen that the limit solution u has a different 
fall off behavior at the origin than the family u{m p ) for 
mP > 0, in particular the conformal factor for m? > 
behaves like 0(l/r) at r — but in the limit m p — > it 
behaves like 0(l/y/r). This implies a change in the fall 
off behavior off the physical fields at the end r = 0. This 
end will not be asymptotically flat in the extreme limit. 

We illustrate the same phenomena with two important 
examples. The first one is the Reissner-Nordstrom black 
hole initial data. In isotropic coordinates, a canonical 
slice t = const for the Reissner-Nordstrom black hole 
spacetime with mass Madm and charge q (for a black 
hole we always have q < Madm) defines the following 
initial data set 



7ij = ip%, Kv=0, r=0, /i=-TT. (14) 

r <p 

where the conformal factor is given by 

<P = ttV(9 + 2r + M AD M)(-q + 2r + Madm)- (15) 
Ir 

The conformal factor satisfies the following equation 
(analog to Eq. {S} ) 



Ap 



Ar 4 ip 3 ' 

For this solution define the parameter 



mP = \/M ADM -q 2 , 
and the function u by 



u = ip — 1 



2r 



(16) 



(17) 



(18) 
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For mP > 0, the function u is bounded at r — 0. The 
extreme limit q — > Madm corresponds to mP — ► 0. In 
this limit we have 



<P 



M A 



DM 



1. 



U = If — 1, 



(19) 



and hence ip = u = 0(\/^/r). Note that although 
Eq. (fT6|) is different from (J5]) , the powers of r and 99 on 
the right hand side are such that if we assume u — 0(r a ) 
at r = we get a = —1/2 as in the Bowen-York spinning 
black hole. 

The second example is given by the Kerr black hole 
initial data in quasi-isotropic coordinates. Let S and 
Madm be the total angular momentum and mass of the 
Kerr spacetimc. Define the parameter vnP by 



mr 



Madm 



S 



M 2 

ADM 



(20) 



As in the previous example, the extreme limit yS — > 
Madm of Kerr metric correspond m p — > 0. The explicit 
expression of these data can be found in [o5j ]. in this 
reference it is proved that the conformal factor behaves 
in a similar way as the Reissner-Nordstrom in the limit 
m p — > 0. For completeness, we reproduce this calcula- 
tion and adapt it to our setting. We use the coordinate 
transformation 66] 



m + a \ I m — a 
r = f 1+ — — 111 



2f 



2r 



(21) 



where r is the standard Boyer-Lindquist radial coordi- 
nate, which puts the Kerr metric in the quasi-isotropic 
form 

ds 2 = <f K [e~ 2<lK (df 2 + f 2 d9 2 ) + f 2 sin B 2 d^ 2 \ , (22) 
where ds 2 is the spatial line element, 



-2q K 



r 2 + a 2 cos 2 9 



„,2 i „2 i 2roa 2 rs in 2 



and 
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2 1 2 

r + a 



r 2 +a 2 cos 2 



2ma 2 r sin 



a 2 cos 2 



The function </> 4 has the expansion 
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} 3 — a 2 m a 2 + 3m 2 



16f 4 2f 3 
2m (-2cos(26»)a 2 + a 5 



(m 2 — a 2 ) r 



2f 2 

+ 0(1) 



when a < m. However, for a — m expansion 
singular. In this case </> has the form: 

4 4m(m + f) (2m 2 + 2fm + f 2 ) 

^ ~ f 2 (cos(26»)m 2 + 3m 2 + 4fm + 2f 2 ) + 1 ' 



(23) 
(24) 

f (25) 
(26) 
25]) is 

(27) 



which has the expansion 

4 8m 2 32m cos 2 (6») 



VK (cos(2(9) + 3)f 2 (cos(26») + 3) 2 f 
Thus the behavior of 4>k changes from 



•0(1). (28) 



X - y/a 2 - m 2 /(2f) + 0(1) 



to 



\Jlmjr 
^l + cos 2 



+ 0(r x / 2 ) 



in the extremal case a — m (the horizon is then mapped 
to the limiting surface f = 0). Hence, just as in the max- 
imal Bowen-York case, the f = coordinate singularity 
corresponds to an end which is not asymptotically flat. 



B. Numerical Test of Highly-Spinning BY Initial 
Data 



We first solve the initial data for one black hole pre- 
sented in section III Al We used a modified version of 
BAM_Elliptic thorn |60l. [67| hi order to solve for the case 
where S/M 2 = 1 (here M denote an scale factor) and 
m p = on a uniform grid with resolution h = 0.0035M 
and outer boundary at 0.64M (here we are only inter- 
ested in the singular behavior, which is not affected by 
inaccurate boundary data). In order to avoid the singu- 
larity itself, we constructed the grid such that the origin 
was located halfway between gridpoints. We fit "f xx along 
the lines y = z = h/2 and x = y = h/2 to the form: 
Ixx = a + b/(x 2 + 2c 2 ) (with a similar form for z) and 
perform a non-linear least-squares fit (See Fig. [T|). We 
find that j xx (x) = 11.483 + 2.5042/(a; 2 + 2(0.00231057) 2 
and j xx (z) = 265.319 + 1.91325/(x 2 + 2(0.00192925) 2 ). 
Note that the c parameter differs from the expected h/2 
but this functional form captures the expected singular 
behavior to high accuracy. In Fig. [5] we plot the function 
u for various choices of mass parameter m p for a config- 
uration consisting of two black holes, one at the origin 
and one at x — — 20M, with the former hole having spin 
S/M 2 = 1 and the latter nonspinning. From the plot one 
can see that, although u is finite, u ends to \/\fr as m p 
tends to zero. 

We are not able to calculate the ADM mass Madm of 
the data, because we needed very high resolution near the 
puncture r = 0, leaving too few point outside the hori- 
zon. Instead, we compute the horizon mass mn given 
by tuh = \J m 2 rl + S 2 / (4m 2 rl .) where m m is the irre- 
ducible mass (it is expected for these kind of data that 
mH < Madm , there exists however no proof of this con- 
jecture). Analogous to the ratio (fTTj) we define the quasi 
local ratio a/ ran, where a = S/uih- For these data 
the maximum possible value of this quantity is given by 
a/m H ~ 0.93 (HSIHi- 

Even the horizon mass run is difficult to resolve at the 
initial surface (the horizon is located at only one point 
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FIG. 1: The "/ xx component of the metric on the initial slice 
for a single spinning, non-boosted black hole located at the 
origin with S = 1M 2 and the puncture mass parameter m p = 
0. Here we plot *y xx versus x (y — z — .00175M) and z 
(x = y = 0.00175M) and fits to j xx = a + b/(x 2 + 2c 2 ) (with 



a similar form for z). Note that y x 
the expected ip oc 1/y/r. 

10 6 r ■ , , , 



oc 1/r is consistent with 



fit (X) 

fit (z) 




FIG. 3: The horizon mass and specific spin for a maximal 
BY black hole. Here S = 1M 2 and m p = 0. Note that the 
spin drops rapidly to about a/mn = 0.934 (and is still drop- 
ping) after the hole absorbs significant mass. The expected 
asymptotic value is a/mn = 0.928. 




FIG. 2: The function u for finite m p in the neighborhood of 
the puncture. Here the data consists of two BH (the second 
BH is located at x/M = —20 and is not shown). The BH 
at the origin has spin parameter S/M 2 = 1. Note that u 
approaches \j\fr as m v tends to zero. 
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r = 0), we need to perform the evolution of the data to 
compute it at later times. 

These data have axial symmetry and then the spin is a 
conserved quantity. One part of the radiation emitted by 
the data will fall into the black hole increasing its area 
and the other part will scape to infinity. Hence the ratio 
a/mn will monotonically decrease during the evolution. 
This is precisely what we observed in the numerical evo- 



lution presented in Fig. [3] 

In order to measure the horizon mass and the quotient 
a/mn in this case, we performed a second unigrid run, 
this time using fisheye [1^, , with a central resolution 
of h = M/64 and outer boundaries at 65M. Initially the 
horizon has a coordinate radius of zero (akin to extreme 
Kerr in quasi-isotropic coordinates) that grows to r ~ 
0.3M at t = 15M as it absorbs the spurious radiation 
produced by the conformally flat initial data. In Fig. [3] we 
show the horizon mass and specific quotient for this run. 
Note the rapid drop off in the ratio a /tuh between t = 
10M and £ = 15. Due to the relatively poor resolution of 
this run, and the short evolution time, it is not clear if 
it asymptotes to a/mn = 0.933, or will continue to drop 
to the numerically predicted value of 0.928. 



C. Spinning-black- hole-binary initial data 

The initial data techniques of Sec. Ill Al can be ex- 
tended to multiple spinning black holes with linear mo- 
mentum [60l ]. Here too the 3-metric on the initial 
slice has the form y a j, — (ipBL + u) 4 S a b, where ipBL 
is the Brill-Lindquist conformal factor, 8 a b is the Eu- 
clidean metric, and u is (at least) C 2 on the punc- 
tures. The Brill-Lindquist conformal factor is given by 
ipBL = 1 + Si=r m ii]/ i^l^ — f[i] I), where the sum is over 
all punctures, mf^ is the mass parameter of puncture i 
(rnP-i is not the horizon mass associated with puncture i), 

and f[i] is the coordinate location of puncture i (we use 
the notation [i] to distinguish the puncture label from the 
tensor indices in the equations below) . The extrinsic cur- 
vature K a b is given by the Bowen-York (BY) ansatz [58| 
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0.5059 



III. NUMERICAL TECHNIQUES 




FIG. 4: The specific spins and horizon masses for a binary 
containing nearly maximal BY spinning holes. The initial 
spins of the two holes are ajmn = 0.967. The individ- 
ual horizon masses increase rapidly near t ~ 15M as the 
black holes absorb the spurious radiation. The spin drops to 
a/ran = 0.924 since the spurious radiation does not increase 
the angular momentum of the holes. 



and has the form K ab = if 10 K ab , where 



j|2( 2P W»W-(^-»M»w)^»Wc) 
^(2S [t]c n lt]d e cd ^ a n b l), (29) 



2|f- f[ 
3 



<P = i^BL+u, ft[i] = (r-r^/lr-r^], and P^ and are 
the linear and angular momenta of puncture i. As shown 
in Sec. Ill Al for a single puncture, if mP = 0, P = 0, 
and S =/= then u is no longer finite at the puncture loca- 
tion, but has a singularity. Under these conditions, 
the resulting black hole will have the maximum possible 
specific spin for Bowen-York type data. 

We compute and evolve the black-hole-binary data de- 
scribed above. The ratio a/ ran is a quasi local quantity 
which is well defined for individual black holes in a binary 
system (in contrast with the ratio £5 in which appears 
the ADM mass which is a global quantity). We can ob- 
tain specific ratios ajmn nearly equal to the maximum 
allowed for BY holes by setting mP sufficiently small. 
For the runs presented below we use the puncture ap- 
proach along with the TwoPUNCTURES [ill thorn 
to compute initial data for black-hole binaries. In Fig. Q] 
we show the isolated horizon spin for the individual holes 
in a near maximally spinning BY binary (configuration 
MR0 described below). Note that the spins drop signifi- 
cantly within 10M - 20M of evolution. 



We evolve the black-hole-binary data-sets using the 
LazEv [9l |7 2fl implementation of the moving puncture 
approach |9i, llOj] . In our version of the moving puncture 
approach we replace the BSSN [73|, El, Ej| conformal 
exponent 0, which has logarithmic singularities at the 
punctures, with the initially C field \ — exp(— 40). This 
new variable, along with the other BSSN variables, will 
remain finite provided that one uses a suitable choice for 
the gauge. An alternative approach uses standard finite 
differencing of (f> (loj . 

We use the Carpet 176] mesh refinement driver to pro- 
vide a 'moving boxes' style mesh refinement. In this ap- 
proach refined grids of fixed size are arranged about the 
coordinate centers of both holes. The Carpet code then 
moves these fine grids about the computational domain 
by following the trajectories of the two black holes. 

We obtain accurate, convergent waveforms and horizon 
parameters by evolving this system in conjunction with 
a modified 1+log lapse and a modified Gamma-driver 
shift condition [9|, [77], and an initial lapse a(t = 0) = 
2/(1 + ip% L ). The lapse and shift are evolved with 



(d t -f3 l di)a = -2aK 
d t (3 a = B a 
d t B a = 3/4c> t f Q 



r)B a , 



(30a) 
(30b) 
(30c) 



These gauge conditions require careful treatment of x, 
the inverse of the three-metric conformal factor, near the 
puncture in order for the system to remain stable d [70l . 
[78j . In Ref. [l2j it was shown that this choice of gauge 
leads to a strongly hyperbolic evolution system provided 
that the shift does not become too large. 

We use AHFinderdirect to locate apparent hori- 
zons. We measure the magnitude of the horizon spin us- 
ing the Isolated Horizon [ll El El [M El El algorithm 
detailed in [80] . This algorithm is based on finding an ap- 
proximate rotational Killing vector (i.e. an approximate 
rotational symmetry) on the horizon, and given this ap- 
proximate Killing vector £ a , the spin magnitude is 



S 



is] 



1 

8^ 



(CR b K ab )d 2 v 



(31) 



AH 



where K a b is the extrinsic curvature of the 3D-slice, d 2 V 
is the natural volume element intrinsic to the horizon, 
and R a is the outward pointing unit vector normal to 
the horizon on the 3D-slice. We measure the direction of 
the spin by finding the coordinate line joining the poles 
of this Killing vector field using the technique introduced 
in [l6| . Our algorithm for finding the poles of the Killing 
vector field has an accuracy of - 2° (see pi] for details). 

We measure radiated energy, linear momentum, and 
angular momentum, in terms of ^4, using the formulae 
provided in Refs. [8l], [82|]. However, rather than using 
the full we decompose it into £ and m modes and 
solve for the radiated linear momentum, dropping terms 
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TABLE I: Initial data parameters for the maximum recoil 
configurations MRO - MR31 5 and the maximum hangup con- 
figuration MH. For MRO - MR3I5 the punctures are lo- 
cated at x/M = (±3.564036838, 0,0), with momenta p/M = 
(0, ±0.1254868859, 0) and spins S = ±(S x ,S y ,Q). In all 
cases the orbital frequency is Muj = 0.045 and the punc- 
ture mass parameters are mP /M — 0.08967. For MH 
the punctures are located at X/M = (±4.083304018,0,0), 
with momenta P = (0, ±0.1046285561, 0) and spin S/M 2 = 
±(0,0,0.2362330001). The orbital frequency is Mlo = 0.035 
and the puncture mass parameters are m p /M = 0.107949. 



Config 


Madm 


S x /M 2 


Sy/M 2 


MRO 


1.000000 





0.23642497 


MR45 


0.999493 


-0.16717770 


0.16717770 


MR90 


0.998982 


-0.23642497 





MR135 


0.999493 


-0.16717770 


-0.16717770 


MR180 


1.000000 





-0.23642497 


MR225 


0.999493 


0.16717770 


-0.16717770 


MR270 


0.998982 


0.23642497 





MR315 


0.999493 


0.16717770 


0.16717770 


MH 


1.000000 







FIG. 5: The xy projections of the puncture trajectories for 
MRO and MR45, the latter rotated by an angle cor so that 
the late inspiral and merger phases overlap. 
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A. Large Recoil Velocities 



with I > 5. The formulae in Refs. [8lJ, |82| are valid 
at r = oo. We obtain highly accurate values for these 
quantities by solving for them on spheres of finite radius 
(typically r/M = 25,30,35,40), fitting the results to a 
polynomial dependence in / = 1 fr, and extrapolating to 
1 = 0. We perform fits based on a linear and quadratic 
dependence on I, and take the final values to be the aver- 
age of these two extrapolations with the differences being 
the extrapolation error. 



IV. RESULTS 

Evolving black holes with specific spins of a /run ~ 
0.92 is challenging because the horizons appear quite 
small. In our coordinates, the initial horizon radii were 
0.04M. We evolved these data using 14 levels of refine- 
ment with a finest resolution of h = Af/320 (the high- 
est resolution reported so far in numerical simulations 
of binary black holes). The outer boundaries were lo- 
cated at ±640M in all directions and the resolution on 
the coarsest grid was h — 12. 8M. We obtained the spin, 
position, and momentum parameters of the initial data 
using third-order post-Newtonian parameters for equal- 
mass quasi-circular binaries with spins aligned with the 
linear momentum of the two holes, and aligned with the 
orbital angular momentum. In all cases we took the spin 
of the two holes to be a/mn = 0.92. We set the punc- 
ture mass parameters of the two holes such that the to- 
tal ADM mass was 1M. The initial data parameters are 
summarized in Table HI Note that we normalize the MRO 
configuration to ADM mass of 1M, but keep the same 
mass parameters when we rotate the spin. Hence MR45 
- MR315 have slightly different ADM masses. 



In Ref. [33j we proposed a semi-empirical formulae for 
the dependence of the merger recoil velocity on the spins 
and mass ratio of the two black holes in a binary. Our 
formula predicts that the largest recoils occur for two 
equal-mass, equal-spin black holes with spins pointing in 
the orbital plane and counter aligned with each other. In 
this configuration the recoil is proportional to the spin- 
amplitude and varies sinusoidally with the angle between 
the direction of the spins at merger and the linear mo- 
mentum direction. The resulting recoil will be perpendic- 
ular to the orbital plane. We were able to test this pre- 
diction (35| by evolving a binary with spins a/mn = 0.5 
(initially pointing in the direction of the linear momen- 
tum) and then evolving a set of binary configuration with 
the same orbital parameter, but with the initial spin di- 
rections rotated by an angle 9. We found that we could 
fit the resulting recoil velocities to a cos(0 — ©o) de- 
pendence to high accuracy. Implicit in this approach is 
that the angle between the spins and linear momentum 
at merger is given by the angle at merger for the = 
configuration plus O. This, in turn, requires that the xy 
projection of the trajectories be independent of 0. 

Unlike in Ref. 35j, here we found for the MRO - 
MR315 configurations that rotating the initial spin di- 
rection alters the ^-projection of the trajectories by in- 
troducing varying ellipticity to the orbit and possibly due 
to spin-orbit coupling. However, after rotating the xy- 
projected trajectories for M45 - M315 by an angle ro t 
(See Fig. and Ref. [83[ for a description of the tech- 
nique), we found that they overlap for the late inspiral 
and plunge with projected MRO trajectory. We then take 
the initial spin orientation, plus this trajectory rotation, 
as the angle between MR45 - MR315 spin direction at 
merger with the MRO spin direction at merger. We then 
fit the ^-component of the recoil versus this angle. The 
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TABLE II: Recoil velocity (in the z-direction) for the MRO — 
MR315 configuration, initial angle O CO n between the spin di- 
rections of the MR45 — MR315 and MRO, and (approximate) 
angle O cor between spin directions of MR45 — MR315 and 
MRO at merger. Here O cor = 6 CO n + Orot- 



Config 


©con 


@cor 


^(kms- 1 ) 


MRO 


0° 


0° 


2372 ± 12 


MR45 


45° 


15° 


2887 ± 27 


MR90 


90° 


40° 


3254 ± 19 


MR135 


135° 


92° 


2226 ± 6 


MR180 


180° 


186° 


-2563 ± 8 


MR225 


225° 


195° 


-2873 ± 23 


MR270 


270° 


205° 


-3193 ±45 


MR315 


315° 


250° 


-2910 ± 1 



FIG. 6: Recoil velocity versus corrected rotation angle for the 
MRO — MR315 configurations and a non-linear least-squares 
fit to a simple sinusoidal behavior. Note that the corrected 
rotation angles are not distributed uniformly in the range 
(0,2tt). 

4000 | . 1 . 1 . 1 . 1 1 




-4000 1 1 1 1 1 1 1 1 L 

12 3 4 
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results are summarized in Table ITT1 and Fig. [6] (with sev- 
eral runs having significantly higher recoil velocities than 
any previous simulation). A fit of v z versus CO r (where 
cor is the corrected angle between MR45 - MR315) and 
MRO gives: !) z (km s" 1 ) = 3290.14 cos(6 cor - 0.765885) 
where O CO r is measured in radians and the confidence in- 
terval for the amplitude is (3243,3337). Our empirical 
formula predicts an amplitude of 3461 ± 58. These re- 
sults are with 2.3 <r of the prediction for a/mn — 0.923. 
A relatively large error is not unexpected due to the dif- 
ficulty in evolving systems with such small scale features 
and the ellipticity introduced by the spurious radiation. 
On the other hand, this may also indicate a non-linear 
(in a i 'mii) term is present that reduces the maximum 
possible recoil. 



FIG. 7: The puncture trajectories and first common appar- 
ent horizon for the MH configuration. Note that the binary 
completes ~ 7.5 orbits prior to merger. 
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B. Orbital Hangup 

Of particular interest is the spin-orbit hangup ef- 
fect [lj, S3, H3 when the two spins are aligned with the 
orbital angular momentum. Here we examine configura- 
tion MH, where the two spins have near maximal (for BY 
data) spins. We evolved configuration MH with 14 levels 
of refinement, maximum resolution of h — Af/320, and 
outer boundaries at 1281A/ . In Fig. [7] we show the xy 
trajectories of the punctures as well as the first common 
apparent horizon. Note that the binary completes ~ 7.5 
orbits before the first common apparent horizon forms. 
In Fig. [5] we show r = \ f\ — fa| (where fi is the location 
of puncture i) versus orbital phase </> rbit- The initial 
eccentricity, as is evident by the oscillation in r(</> or bit), 
damp with time. The hangup effect is clearly seen in 
Fig. [5J which shows the xy projections of the trajectory 
difference f\ — f*2 for the MH and MRO configurations. 
In Fig. [10] we show the real part of the (£ — 2, m = 2) 
mode of ip4 for the MH configuration. We measure the 
remnant mass and spin three different ways: from the 
isolated horizon formalism, from the radiated energy and 
angular momentum, and from the quasi-normal decay of 
the late-time waveform. Results for the isolated horizon 
calculation were affected by late-time boundary effects 
due to relatively poor resolution in the outer zones. The 
results are summarized in Table Mil To calculate the 
horizon mass and spin from the quasi-normal modes we 
used the results of [85| and a fit to exp(— t/r) sin(a>£ — 9q) 
for the real and imaginary parts of the (£ = 2, m = 2) 
mode of ^4 at r = 40M. We found t/M = 15.097±.0074 
and Muj = 0.76579 ±0.00006 (the errors reported are the 
differences in r and u> for the real and imaginary parts of 
tpi). We also calculated the remnant mass and spin from 
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FIG. 8: The puncture separation, and first derivative (see 
inset), versus orbital phase for the MH configuration. Note 
the decaying oscillations that indicate the ellipticity is signif- 
icantly reduced after 2 orbits. 




FIG. 10: The real part of the (( = 2, m = 2) mode of ip4, 
for the MH configuration showing the orbital dynamics and 
quasi-normal decay. 
0.2 



3 4 5 

<P„ rbi /(2 it) 




1000 



FIG. 9: The xy projection of the trajectory difference ri — r*2 
for the MH and MR0 configurations. Note the significant 
hangup effect when the spins are parallel to the orbital angu- 
lar momentum (MH). 




-8-3 2 7 



(x, - x 2 )/M 



the radiated energy and angular momentum. Here too, 
inaccuracies due to relatively poor resolution in the outer 
zones affect the calculation. For all three methods, the 
final spin is in qualitative agreement with the prediction 
of a/rriH — 0.928 in but differs from the predic- 

tion a/m H = 0.9400 ± 0.0019 in [H for initial specific 
spins of 0.92. 



TABLE III: Remnant horizon mass and spin for configuration 
MH based on the isolated horizon (IH) calculation, the radi- 
ated energy and angular momentum, and the quasi-normal 
frequency (QNM). 



IH 


Radiation 


QNM 


m H 0.9095 ± 0.0005 
a/m H 0.922 ±0.001 


0.9162 ± 0.0024 
0.928 ±0.015 


0.9146 ±0.0002 
0.922 ±0.001 



V. CONCLUSION 

We have evolved equal-mass, equal-spin black-hole- 
binary configuration with nearly maximal BY spin of 
a/m,H = 0.923, the highest spins simulated thus far, 
both for spins pointing in the same direction as the or- 
bital angular momentum and for spins (counter-aligned) 
pointing in the orbital plane. In the former case we 
see a significant orbital hangup (for seven orbits prior 
to merger), and confirmed that the remnant spin is 
non-maximal and agrees with our previous predictions 
based on a least-squares fit of remnant spin versus initial 
spin 0, HH, Gil HH • While in latter case we find that the 
maximum recoil is iWoii = 3290 ± 47 km s , in qualita- 
tive agreement with our prediction of 3461 ± 58 km s _1 . 
The deviation of our measured recoil velocity from the 
predicted value is likely due to the relatively poor ef- 
fective resolution (i.e. the number of gridpoints across 
the initial horizons), as well as eccentricities introduced 
by the significant amount of spurious radiation or non- 
linear corrections to the kicks formula. We have con- 
firmed the sinusoidal dependence of the recoil on the ini- 
tial spin direction and that the recoil varies essentially 
linearly with the magnitude of the spin for a fixed ini- 
tial spin direction. Note that the measured recoil of 
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"^recoil = 3254 ± 19 km s _1 for the MR90 configuration 
is the largest recoil velocity measured for an actual sim- 
ulation to-date, and would expel the BH remnant from 
any known galaxy. 

All our numerical calculations confirm cosmic censor- 
ship. In particular, highly spinning black holes, which are 
close to the extreme Kerr limit, do not develop naked sin- 
gularities. Moreover, the system always decays asymp- 
totically in time to a final black hole which satisfy the 
inequality between mass and spin of the Kerr black hole. 

There is a curious analogous behavior between the con- 
formal factor of the extreme spinning BY initial data as 
r — * 0, i.e. tp ~ and the late-time behavior of the 

determinant of the 3-metric for a Schwarzschild black 
hole [13, [H, H^, [9(| with the standard moving-puncture 
choice for the gauge Eqs. (f3"0"]) . Here r — corresponds 
to the horizon on the initial slice, which has finite sur- 
face area, while in the case of Schwarzschild, if we take 
the limit t — > oo and then r — ► 0, we approach a sphere 
of finite surface area inside the horizon. It seems that, 
as we increase the spin of the black hole to its maxi- 
mum allowed value, we find a new, extreme solution that 
has a different behavior from the non-extreme cases (1/r 



vs. l/^fr). The limit of stationarity, for t — > oo of sub- 
maximal data, also leads to a new behavior for the con- 
formal factor <j> ~ 1/v^j n °t present at finite time. In 
both cases the data transition from a slicing that con- 
tains two asymptotically flat ends to one that contains 
one asymptotically flat end and one cylindrical end. 
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